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We study the nature of the ground state of the frustrated Ji — J2 model and the Ji — J3 model 
using a variational algorithm based on projected entangled-pair states (PEPS). By investigating 
spin-spin correlation functions, we observe a separation in regions with long-range and short-range 
order. A direct comparison with exact diagonalizations in the subspace of short-range valence bond 
singlets reveals that the system is well described by states within this subset in the short-range 
order regions. We discuss the question whether the system forms a spin-liquid, a plaquette valence 
bond crystal or a columnar dimer crystal in these regions. 

PACS numbers: 75.10.-b, 75.10.Jm, 75.40.Mg, 02.70.-c, 03.67.-a 



I. INTRODUCTION 

Frustrated spin-systems have attracted a lot of inter- 
est in the last years, because they may possess exotic 
ground states that are very different from conventional 
Neel-ordered states. Such states are especially intrigu- 
ing, as connections to high-T^ superconductivity have 
been put forward^. They are usually characterized by a 
break down of long-range order: the system reorganizes 
in a quantum state where only local antiferromagnetic 
correlations are present. The class of such states, named 
Short Range Valence Bond States (SRVB), encompasses 
a broad range of phases: they range from valence bond 
crystals (VBC) with broken translational symmetry to 
pure spin liquids that have all symmetries restored. 

Studies of frustrated systems are especially challeng- 
ing, because Quantum Monte-Carlo (QMC) studies 
are hindered by the sign-problem-^ and Density Matrix 
Renormalization Group (DMRG)^i^i^ investigations are 
hard to perform for systems with dimensionality larger 
than one. Other methods that have been developed 
to take on these systems are, for example, the coupled 
cluster method^, DMRG combined with QMC- and ex- 
act diagonalizations within the subspace of SRVB^*^. In 
this paper, we would like to give the recently developed 
PEPS-algorithm^°'^^ a try. This algorithm has already 
been successfully applied to the Heisenberg antiferro- 
magneljiSi, the Bose-Hubbard modelii and the frustrated 
Shastry-Sutherland modeli^. 

In the following, we focus on the J1 — J2 — J3 model. We 
discuss our observations and possible implications. Be- 
fore that, we give a quick review of the PEPS algorithm in 
section im For a thorough explaination of the algorithm, 
we refer to Refs. ilftJlulS . With the implementation of 
the algorithm, we follow in large part Ref. [ll|. 

II. THE PEPS-ALGORITHM REVISITED 

The PEPS algorithm is a variational algorithm within 
the class of Projected Engangled Pair States (PEPS) ^"'^^ . 
These states are expedient for a numerical study of 
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FIG. 1: Classical phase diagram of the Ji — J2 — Js model 
consiting of (I) a Neel phase, (II) a phase consisting of two 
independent sublattices, (III) a helical phase with order {q, q) 
and (IV) a hehcal phase with order (g, tt). The dashed line 
denotes the boundary of the Neel phase when quantum fluc- 
tuations are taken into account. 



many-particle systems: on the one hand, they fulfill the 
area-law that is fundamental to non-critical systems. On 
the other hand, they possess a number of parameters that 
scales polynomially with the system size. This feature 
makes numerical simulations feasible. 

The structure of PEPS is inspired from the structure 
of Matrix Product States (MPS)^'''^^ that form the basis 
for the successful DMRG algorithm^iii^: to each lattice 
site, a tensor is associated that possesses a physical index 
and a certain number of virtual indices. The virtual in- 
dices are contracted according to a scheme that mimics 
the underlying lattice-structure. The dimension of the 
physical index equals the physical dimension of the par- 
ticle residing on that lattice site. The dimension of the 
virtual indices, D, is the internal refinement parameter 
of PEPS: D = 1 specializes the PEPS to a product state; 
the choice D — c?^ (with N being the total number of 
lattice sites and d the physical dimension of one particle) 
enlarges the space of PEPS to the complete Hilbert-space 
of the system. 
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The main idea of the algorithm is to optimize the ten- 
sors such that the energy tends to a minimum. This can 
be done by sequentially optimizing the tensors or by cool- 
ing the system with the simulation of an imaginary time 
evolution. Here, wc follow the latter. Time and memory 
required for this method turn out to be polynomial both 
in TV and D: time scales as N^D^^ and memory as ND^. 

This polynomal scaling allows us to investigate systems 
of sizes that are out of reach with other algorithms. We 
discuss our observations in the following sections. 



III. J1-J2- J3 MODEL 

In the J1 — J2 — J3 model on a square lattice, frustration 
is caused by the competition between first, second and 
third neighbor interactions of magnitudes Ji, J2 and J3 
respectively: 



The phase diagram of this model is involved and still 
controversal. Of special interest are the regimes of maxi- 
mal frustration that are suspected of having non-classical 
ground states. 

Let us first review the classical limi t^^i^^i^^i^^ (5* — > 
00). In this limit, the system possesses four phases, 
as shown in figure [TJ the usual Neel phase, two spi- 
ral antiferromagnetic phases ordered at {q, q) and (g, tt) 
and a phase consisting of two independent sublattices. 
The Neel phase is bounded by the classical critical line 
( J2 + 2 J3)/ Ji = 1/2. 

When quantum fluctuations are taken into account, 
the phase diagram changes considerabl y^i^^i^°i^^i^^ : the 
Neel phase substantially extends to larger values of J3, 
up to the line of maximal frustration ( J2 + J3)/ Ji — 1/2. 
In the vicinity of this line, it is believed that the classical 
ordered ground state is destabilized and a singlet ground 
state is formed^. The precise nature of this state is still 
controversal: suggestions include columnar valence bond 
crystals^'^, plaquette states.* and spin hquid a^'^i^^'^^i^^'^* . 
Special attention has been devoted to the end-points of 
the line: the point at J3/J1 = 1/2 that separates the 
Neel from the spiral antiferromagnetic phase^^^i^^, and 
the tri-critial point at J2/J1 = 1/2 at which 3 phases 
meet2>25^>S9. 

In the following, we focus on two lines in the phase 
diagram that include these points: J2 = (Ji — J3 
model) and J3 = (Ji — J2 model). In both cases, we 
apply the PEPS-algorithm and discuss our observations. 
Our investigations are thereby inspired by the discussions 
in Ref. i8|. We obtain the PEPS-approximation of the 
ground state via imaginary time evolution, as described 
previously. In the course of the time evolution, we in- 
crease the virtual dimension D, starting from D = 2, and 
decrease the Trotter-step St, starting from St — 0.03/ Ji, 
until the relative changes in the ground state energy are 
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FIG. 2: Energy per site as a function of the imaginary time 
evolution steps. The lattice size is 10 x 10 and J3/J1 = 0.5. 
The inset shows the magnified tail of the energy function. 
The vertical lines separate the plot into three regions: in these 
regions, the Trotter-step St is chosen equal to 0.03/ Ji , 0.01 / Ji 
and 0.003/Ji respectively. 



of order 10^'^. This is usually achieved for 13 = 4 and 
St — 0.003/Ji. Typically, D = 4 only leads to a minor 
reduction of the energy. The energy as a function of the 
imaginary time evolution steps is plotted in figure [5] for 
the special case of a 10 x 10 lattice and J3/J1 — 0.5. 



A. Spin— Spin Correlations 

In order to get a first idea of the nature of the ground 
state, we calculate the spin-spin correlation functions 
{sk ■ si) and the corresponding static structure factor. 



(1) 



kl 



in the PEPS-approximation of the ground state for < 
J3/J1 < 1 at J2 = and < J2/J1 < 1 at J3 = 0. 



Ji - J3 Model 

In the case of J2 = 0, we observe that S'(q) is peaked 
at q = (tt, tt) for J3/J1 < 0.5, indicating long-range Neel 
order. This order disappears at J3/ Ji 0.5 at which the 
structure factor becomes smooth. At around J3/J1 ~ 0.7 
peaks at (±7r/2, ±7r/2) reappear, indicating a revival of 
incommensurate long-range order. Figure [3] depicts the 
functional characteristics of the structure factor S'(q) at 
selected points J3/J1, calculated with a, D = 3 PEPS on 
a 14 X 14 lattice. 

The collapse of long-range order is confirmed by a 
direct observation of the spin-spin correlations. These 
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FIG. 3: Structure factor S'(q) for J3/J1 = 0, 0.2, 0.5 and 0.9. 
The results were obtained for a 14 x 14-lattice and virtual 
dimension D — 3. 



FIG. 5: Structure factors S{'k, tt) (upper plot) and 
S'(7r/2, 7r/2) (lower plot) as functions of J3/J1 for various 
particle-numbers A'^. The solid lines represent extrapolations 
to the thermodynamic limit (see text). 
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FIG. 4: Spin-spin correlations {si ■ Sj) as functions of the 
distance A — for the 14x 14 Ji — J3 model at J3/J1 = 0, 

0.2, 0.5 and 0.9. The results have been obtained using a 
D = 3-PEPS Ansatz. 

are shown in figure [H Here, correlations with the cen- 
tral spin are plotted as functions of the distance. As 
can be seen, the spins are antiferromagnetically ordered 
for J3/J1 < 0.5. For J3/J1 > 0.5, every second spin 
possesses antiferromagnetic order. However, the long- 
range order of the spins disappears in the vicinity of 
J3/J1 - 0.5. 

The separation in long-range and short-range or- 
der regimes gets more and more evident with growing 
particle-number. This can be gathered from figure [5] 
Here, 5'(7r,7r) and S{tt/2,tt/2) are plotted as functions of 
J^l Jii evaluated for various particle-numbers N . An ex- 
trapolation to the thermodynamic limit {N 00) shows 



FIG. 6: Structure factors S{n, n) (upper plots) and 
S(7r/2,7r/2) (lower plots) as functions of 1/A'^ for selected 
points Ja/ Ji . The dotted lines represent the best fits by poly- 
nomials of degree 3 in 1/A*'. 



that S{'k,7t) remains finite within the region J3/J1 < 0.3 
and is zero otherwise. On the other hand, 5(7r/2,7r/2) 
is zero up to J3/J1 ~ 0.7 and finite for larger values of 
J1/J3. The region of short-range order is thus narrowed 
down to 0.3 < J3/ Ji < 0.6. 

The extrapolation was obtained by fitting a polynomial 
of degree 3 in 1/7V to values of S'(q) for various values 
of N. The error bars in the plot indicate the error es- 
timates for the predictions. The data points are usually 
fitted very well with such a polynomial - as can be seen 
in figure [HI In this figure, for selected points J3/J1, the 
scaling of S{q) with N is plotted and the best fit by a 
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FIG. 7; Structure factor S'(q) for J2/J1 equal to 0, 0.3, 0.5 
and 0.9. The results were obtained for a 14 x 14-lattice and 
D = 3. 

polynomial of 3rd degree in 1/N is drawn. From this fit, 
we obtain the predictions for the thermodynamic hmit. 



Ji - J2 Model 

For J3 = 0, the static structure factor shows the fol- 
lowing behavior: for J2/J1 < 0.5 the structure factor 
is peaked at (tTjTt) which indicates long-range Neel or- 
der. For J2/J1 larger than 0.5, columnar long-range or- 
der develops which is detected by a peak at (0, tt). In 
fact, this columnar long-range order reveals an order- 
by-disorder phenomenon^: quantum fluctuations select 
from the huge manifold of classical ground states configu- 
rations where all spins are parallel in a given direction. In 
the regime in-between the Neel and the columnar phase 
the peaks disappear and long-range order breaks down. 
The structure factor obtained from the PEPS-calculation 
for D = 3 and a 14 x 14 lattice is plotted in figure [7] for 
selected values of J2 / Ji . 

In figure[8l S{tt, tt) and S{Tr, 0) are plotted as functions 
of J2/J1 for various particle-numbers N. The extrapo- 
lation to the thermodynamic limit indicates that S{n, tt) 
remains finite for J2/J1 < 0.5, and that S'(7r, 0) is fi- 
nite for J2/J1 > 0.5. Thus, the systems consists of two 
phases with long-range order plus a possibly very small 
short-range order phase in the vicinity of J2/J1 ^ 0.5. 

IV. SHORT-RANGE RESONATING VALENCE 
BOND STATES 

The considerable decrease of the correlation length at 
J3/J1 ~ 0.5 and J2/J1 ^ 0.5 opens the possibihty for 
a short-range resonating valence bond state (SRVB) in 
this area. We investigate this possibility by doing a direct 
comparison of the PEPS results to results obtained by an 
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FIG. 8; Structure factors S{tv, tt) (upper plot) and S{n, 0) 
(lower plot) as functions of J2/J1 for various particle- 
numbers N. The solid lines represent extrapolations to the 
thermodynamic limit. 



exact diagonalization of the Hamiltonian in the subspace 
of SRVB. 



Ji - J3 Model 

In case of the Ji — J3 model, the overlap between the 
SRVB and the PEPS with virtual dimensions D = 3 on 
a 6 X 6-lattice can be gathered from frgureO (lower plot). 
As it can be seen, the overlap increases up to 99% at the 
point J3/J1 = 0.5 and is significantly smaller in other 
regions. A comparison of energies, however, reveals that 
the set of valence bond states does not cover all terms in 
the ground state. As shown in the upper plot of figure [51 
the energies of the diagonalization within the subspace 
of SRVB are - though very close to the PEPS-results at 
J3/J1 ^ 0.5 - always higher than the energies obtained 
within the set of PEPS. Thus, the true ground state at 
J3/J1 ~ 0.5 might contain a small fraction of valence 
bond terms that have longer range. 

Even if it is assumed that the ground state is a pure 
SRVB in the vicinity of J3/ Ji — 0.5, its properties can be 
very rich and it needs a more precise classification. On 
the one hand, it could be a state with broken transla- 
tional symmetry - such as a columnar valence bond crys- 
tal or a plaquette state. On the other hand, an equally 
weighted superposition of valence bond states with re- 
stored translational symmetry known as spin-liquid is 
possible. 

The observation of nearest-neighbor spin-spin correla- 
tions {si-Sj) gives us an indication of a plaquette state^i. 
These correlations are shown in figure \W\ for a 8 x 8 lat- 
tice and D — 3. In case of a pure plaquette state, the 
nearest-neighbor spin-spin correlations would be equal 
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FIG. 9: upper plot: Ground state energy of the Ji — J3 
model on a 6 X 6 lattice as a function of J3/J1, obtained by 
diagonalizing within the SRVB subspace (solid line) and by 
PEPS calculations with D = 3 (green dots) and D = A (red 
dots), lower plot: Overlap between the SRVB ground state 
and the D = 3-PEPS ground state. 



FIG. 11: upper plot: Ground state energy of the Ji — J2 
model on a 6 X 6 lattice as a function of J2/J1, obtained by 
diagonalizing within the SKVB subspace (solid line) and by 
PEPS calculations with D = 3 (green dots) and D = 4 (red 
dots), lower plot: Overlap between the SRVB ground state 
and the D = 3-PEPS ground state. 
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FIG. 10: Nearest-neighbor spin-spin correlations at J3/J1 = 
0.5, calculated with D = 3-PEPS. 



to —1/2 on a plaquette and between two plaquettes. In 
our case, the values of the spin-spin correlations deviate 
slightly from these values, nontheless a clear plaquette 
structure remains visible. 

The state we observe has obviously a broken trans- 
lational symmetry. The plaquettes are formed between 
sites (i, j), (j + 1, j), + 1), (i + 1, j + 1) with i and j 
always being odd. The reason for this symmetry break- 
ing are the chosen open boundary conditions and the 
even number of sites in x- and y-direction. The system 
chooses the configuration with a maximal number of pla- 
quettes. In case of open boundary conditions and an even 
number of sites in each direction, this corresponds to the 
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FIG. 12: Nearest-neighbor spin-spin correlations at J2/J1 
0.6, calculated with D = 3-PEPS. 



configuration of plaquettes with i and j being odd. 



Ji - J2 Model 

The vicinity of the ground state to the class of SRVB 
in the region J2/J1 ^ 0.5 is analyzed in figure [TT] The 
lower plot shows the overlap of the PEPS ground state 
with the ground state obtained by exact diagonalization 
within the PEPS subspace. Lattice size 6x6 and virtual 
dimensions D = 3 were considered. The overlap clearly 
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FIG. 13: Plaquette order parameter Qap-yS at J3/J1 = 0.5, 
evaluated on a 8 x 8 lattice. The used virtual dimension is 
D = 3. 



reaches the maximum at J2I J\ — 0.5 and assumes a value 
of about 70%. A comparison of the energies, as shown in 
the upper plot of figure [TI] uncovers that the true ground 
state will not be exactly in the subspace of SRVB: the 
energies of the PEPS calculations are slightly lower than 
the ones obtained from SRVB, even at the critical point. 
Nonetheless, the distance to the subspace of SRVB might 
be very small for J2/J1 ^ 0.5. 

A more precise classification of the ground state turns 
out to be more difficult than in the case of the J\ — J3 
model: observations of the nearest-neighbor spin-spin 
correlations (s^ • Sj ) yield faint indications for a columnar 
dimer state^^. The magnitudes of the nearest-neighbor 
spin-spin correlations for a 8 x 8 lattice, D — i and 
J2/J1 = 0.6 are shown in figure [T21 In a pure dimer 
state, the spin-spin correlations are equal to —3/4 on a 
dimer bond and zero between two dimer bonds. Even 
though these values are not attained, a columnar order 
of dimers is visible. 



FIG. 14: Structure factors Scoi and Svbc as functions of 
J3/J1 for lattice sizes 6x6, 8x8 and 10 x 10, calculated using 
PEPS with D = 3. The solid lines represent extrapolations 
to the thermodynamic limit. 
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FIG. 15: Structure factors Scot and Svbc as functions of 
J2/J1 for lattice sizes 6x6, 8x8 and 10 x 10, calculated using 
PEPS with D = 3. The solid lines represent extrapolations 
to the thermodynamic limit. 



V. OTHER ORDER PARAMETERS 

In order to shed more light on the properties of the 
ground state in the maximally frustrated regime, more 
elaborate order parameters are investigated like the pla- 
quette order parameter—, the columnar order parameter 
or the VBC order parameter—. 

The plaquette order parameter— distinguishes clearly 
a Neel ordered phase from a plaquette valence bond crys- 
tal. It is defined, using the cyclic permutation operator 
Pap-yS of the four spins a, (3, 7 and 5 on one plaquette. 



as 

= 2 (Sq • S/3S^ • S5 + Sq ■ S^S^ ■ S~f - Sa ■ SjSp ■ Sg) 

+1/2 ( 

-1-1/2 (sa ■s^ + Sf)-ss + 1/4) . 

In case of a pure plaquette state, this order parameter 
assumes the value 1 on each plaquette; between the pla- 
quettes, its expectation value is 1/8. The order param- 
eter vanishes in case of the Neel state which lacks the 
cyclic permutation symmetry. 
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The columnar order parameter and the VBC order pa- 
rameter— are defined as structure factors of dimcr-dimer 
correlations: 

Here, A: and I run over all pairs of neighboring sites, i 
and j have fixed values and they point to two neighbor- 
ing sites at the center of the lattice. A represents either 
"col" or "VBC" and the phase factor e\{k,l) is defined 
according to figure 7 in Ref. Nb denotes the num- 
ber of terms with e\{k,l) ^ 0. As discussed in Ref. [Sj, 
Sy BC and Scoi show a different behavior in the plaquette 
and the columnar dimer phase: in the plaquette phase, 
Sv BC remains finite in the thermodynamic limit and Scoi 
tends to zero, whereas in the columnar dimer phase, both 
SvBC and Scoi remain finite. 

All these order parameters are based on dimer-dimer 
correlations which are difficult to obtain with high pre- 
cision from PEPS-calculations. Nonetheless, our results 
give a qualitative picture of the ground state. 



and gives strong indications for a plaquette ordered state 
in the regime of strong frustration. In this regime, our 
results are consistent with those of Ref. d. 



Ji - J2 Model 

In figure [TSl the structure factors Svbc and Scoi as 
functions of J2/J1 are plotted for various lattice sizes. 
The virtual dimension used for the PEPS-simulation was 
D — i. Again, the solid lines represent extrapolations to 
the thermodynamic limit. As can be seen, both quanti- 
ties feature regions in which they are finite and regions 
in which they tend to zero rapidly. In the case of Scoii 
the structure factor is very small for J2/J1 < 0.5 and 
finite otherwise. On the other hand, Svbc is finite for 
J1I J\ ^ 0.7 and close to zero otherwise. This indicates 
the existence of a columnar dimer phase within the region 
0.5 < J2/ Ji < 0.7. 



Ji - J3 Model 

In figure [131 we plot the plaquette order parameter 
evaluated on all square clusters for J3/J1 = 0.5 and a 
8x8 lattice. The results are obtained using a PEPS- 
Ansatz with D = 3. A clear plaquette structure is visi- 
ble, though the absolute values of the expectation values 
differ slightly from the optimal ones. 

In figure [T31 both structure factors Svbc and 

Scoi are 

plotted as functions of J3/J1 for various lattice sizes. 
The values were obtained from a PEPS-simulation with 
D — i. The solid lines indicate extrapolations to the 
thermodynamic limit {N —> 00). The extrapolations 
were obtained by fitting polynomials of 1st degree in 
1/iV. It can be seen that Svbc is clearly peaked in the 
region 0.4 < J1/J3 ^ 0.8. In addition, Scoi is 3 orders 
of magnitude smaller than Svbc in this region. This 
gives us an indication for a plaquette state in the region 
0.4 < J1/J3 < 0.8. 

Thus, the PEPS algorithm reproduces well the proper- 
ties of the Ji — J3 model in the regimes of weak frustration 



VI. CONCLUSIONS 

In conclusion, we have applied the PEPS algorithm to 
the Ji — J2 and the Ji — J3 model and discussed our 
observations. In both models we observed a separation 
in long-range order and short-range order regions. In 
the short-range order regions, the ground state seems to 
lie within the subspace of SRVB-states. Due to lack of 
precision we cannot settle the question of whether the 
ground state is a VBC or a spin-liquid, but there are 
strong indications that the ground state is a plaquette 
state in case of the Ji — J3 model and a columnar dimer 
state in case of the Ji — J2 model. Simulations with 
higher D and the inclusion of symmetries might lead to 
a more concrete answer. 
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